log using "$Rep_smokelabor/2_analysis/output_logs/figure3.log", replace 

	use "$Rep_smokelabor/1_build/regdata/county_quarter.dta", clear 
	tsset fe_countyqtroy rfrnc_yr
	
	** panel a 
	if 1 {
		* static 
		reghdfe pm25 hms_deep [aw=seer_pop], a(fe_countyqtroy fe_styr) vce(cluster countyfip fe_stqtros)
		local s_b=_b[hms_deep]
		local s_l=_b[hms_deep] - invttail(e(df_r),0.025)*_se[hms_deep]
		local s_u=_b[hms_deep] + invttail(e(df_r),0.025)*_se[hms_deep]
		
		* dynamic 
		reghdfe pm25 F2.hms_deep F1.hms_deep hms_deep L1.hms_deep L2.hms_deep [aw=seer_pop], a(fe_countyqtroy fe_styr) vce(cluster countyfip fe_stqtros)
		local d_b_m2=_b[F2.hms_deep]
		local d_l_m2=_b[F2.hms_deep] - invttail(e(df_r),0.025)*_se[F2.hms_deep]
		local d_u_m2=_b[F2.hms_deep] + invttail(e(df_r),0.025)*_se[F2.hms_deep]
		
		local d_b_m1=_b[F.hms_deep]
		local d_l_m1=_b[F.hms_deep] - invttail(e(df_r),0.025)*_se[F.hms_deep]
		local d_u_m1=_b[F.hms_deep] + invttail(e(df_r),0.025)*_se[F.hms_deep]
		
		local d_b_p0=_b[hms_deep]
		local d_l_p0=_b[hms_deep] - invttail(e(df_r),0.025)*_se[hms_deep]
		local d_u_p0=_b[hms_deep] + invttail(e(df_r),0.025)*_se[hms_deep]
		
		local d_b_p1=_b[L.hms_deep]
		local d_l_p1=_b[L.hms_deep] - invttail(e(df_r),0.025)*_se[L.hms_deep]
		local d_u_p1=_b[L.hms_deep] + invttail(e(df_r),0.025)*_se[L.hms_deep]
		
		local d_b_p2=_b[L2.hms_deep]
		local d_l_p2=_b[L2.hms_deep] - invttail(e(df_r),0.025)*_se[L2.hms_deep]
		local d_u_p2=_b[L2.hms_deep] + invttail(e(df_r),0.025)*_se[L2.hms_deep]
		
		* plot 
		preserve 
			clear 
			set obs 5
			
			gen eqtr=_n-3
			foreach var in d_b d_l d_u {
				gen `var'=. 
				replace `var' = ``var'_m2' if eqtr==-2 
				replace `var' = ``var'_m1' if eqtr==-1 
				replace `var' = ``var'_p0' if eqtr== 0 
				replace `var' = ``var'_p1' if eqtr== 1 
				replace `var' = ``var'_p2' if eqtr== 2 
			}
			
			foreach var in s_b s_l s_u {
				gen `var'=. 
				replace `var'=``var'' if eqtr==0
			}
			
			tw 	connected d_b eqtr, col(teal) lw(0.5) msymbol(circle_hollow) || ///
				line d_l d_u eqtr, col(gs9 gs9) lp(dash dash) || ///
				rcap s_l s_u eqtr, col(blue) || ///
				scatter s_b eqtr, col(blue) ///
				xtitle("Event year",size(huge)) ytitle("PM2.5",size(huge)) ///
				ylab(,labsize(vlarge) nogrid) yline(0,lp(dot) lw(0.6) lcol(black)) ///
				xlab(,labsize(vlarge)) ///
				legend(ring(0) pos(2) col(1) size(large) order(1 "dynamic" 5 "static")) ///
				graphregion(color(white)) xsize(4.5)
			gr export "$Rep_smokelabor/2_analysis/output_figures/figure3_a.pdf", replace 
		restore 
	
	}
	
	** panel b 
	if 1 {
		
		* static 
		reghdfe d_pc_qwi_payroll hms_deep [aw=seer_pop], a(fe_countyqtroy fe_styr) vce(cluster countyfip fe_stqtros)
		local s_b=_b[hms_deep]
		local s_l=_b[hms_deep] - invttail(e(df_r),0.025)*_se[hms_deep]
		local s_u=_b[hms_deep] + invttail(e(df_r),0.025)*_se[hms_deep]
		
		* dynamic 
		reghdfe d_pc_qwi_payroll F2.hms_deep F1.hms_deep hms_deep L1.hms_deep L2.hms_deep [aw=seer_pop], a(fe_countyqtroy fe_styr) vce(cluster countyfip fe_stqtros)
		local d_b_m2=_b[F2.hms_deep]
		local d_l_m2=_b[F2.hms_deep] - invttail(e(df_r),0.025)*_se[F2.hms_deep]
		local d_u_m2=_b[F2.hms_deep] + invttail(e(df_r),0.025)*_se[F2.hms_deep]
		
		local d_b_m1=_b[F.hms_deep]
		local d_l_m1=_b[F.hms_deep] - invttail(e(df_r),0.025)*_se[F.hms_deep]
		local d_u_m1=_b[F.hms_deep] + invttail(e(df_r),0.025)*_se[F.hms_deep]
		
		local d_b_p0=_b[hms_deep]
		local d_l_p0=_b[hms_deep] - invttail(e(df_r),0.025)*_se[hms_deep]
		local d_u_p0=_b[hms_deep] + invttail(e(df_r),0.025)*_se[hms_deep]
		
		local d_b_p1=_b[L.hms_deep]
		local d_l_p1=_b[L.hms_deep] - invttail(e(df_r),0.025)*_se[L.hms_deep]
		local d_u_p1=_b[L.hms_deep] + invttail(e(df_r),0.025)*_se[L.hms_deep]
		
		local d_b_p2=_b[L2.hms_deep]
		local d_l_p2=_b[L2.hms_deep] - invttail(e(df_r),0.025)*_se[L2.hms_deep]
		local d_u_p2=_b[L2.hms_deep] + invttail(e(df_r),0.025)*_se[L2.hms_deep]
		
		* plot 
		preserve 
			clear 
			set obs 5
			
			gen eqtr=_n-3
			foreach var in d_b d_l d_u {
				gen `var'=. 
				replace `var' = ``var'_m2' if eqtr==-2 
				replace `var' = ``var'_m1' if eqtr==-1 
				replace `var' = ``var'_p0' if eqtr== 0 
				replace `var' = ``var'_p1' if eqtr== 1 
				replace `var' = ``var'_p2' if eqtr== 2 
			}
			
			foreach var in s_b s_l s_u {
				gen `var'=. 
				replace `var'=``var'' if eqtr==0
			}
			
			tw 	connected d_b eqtr, col(teal) lw(0.5) msymbol(circle_hollow) || ///
				line d_l d_u eqtr, col(gs9 gs9) lp(dash dash) || ///
				rcap s_l s_u eqtr, col(blue) || ///
				scatter s_b eqtr, col(blue) ///
				xtitle("Event year",size(huge)) ytitle("Income per capita",size(huge)) ///
				ylab(,labsize(vlarge) nogrid) yline(0,lp(dot) lw(0.6) lcol(black)) ///
				xlab(,labsize(vlarge)) ///
				legend(ring(0) pos(5) col(1) size(large) order(1 "dynamic" 5 "static")) ///
				graphregion(color(white)) xsize(4.5)
			gr export "$Rep_smokelabor/2_analysis/output_figures/figure3_b.pdf", replace 
		restore 
	
	}
			
	** panel c 
	if 1 {
		
		* static 
		reghdfe d_pmil_qwi_emptotal hms_deep [aw=seer_pop16plus], a(fe_countyqtroy fe_styr) vce(cluster countyfip fe_stqtros)
		local s_b=_b[hms_deep]
		local s_l=_b[hms_deep] - invttail(e(df_r),0.025)*_se[hms_deep]
		local s_u=_b[hms_deep] + invttail(e(df_r),0.025)*_se[hms_deep]
		
		* dynamic 
		reghdfe d_pmil_qwi_emptotal F2.hms_deep F1.hms_deep hms_deep L1.hms_deep L2.hms_deep [aw=seer_pop16plus], a(fe_countyqtroy fe_styr) vce(cluster countyfip fe_stqtros)
		local d_b_m2=_b[F2.hms_deep]
		local d_l_m2=_b[F2.hms_deep] - invttail(e(df_r),0.025)*_se[F2.hms_deep]
		local d_u_m2=_b[F2.hms_deep] + invttail(e(df_r),0.025)*_se[F2.hms_deep]
		
		local d_b_m1=_b[F.hms_deep]
		local d_l_m1=_b[F.hms_deep] - invttail(e(df_r),0.025)*_se[F.hms_deep]
		local d_u_m1=_b[F.hms_deep] + invttail(e(df_r),0.025)*_se[F.hms_deep]
		
		local d_b_p0=_b[hms_deep]
		local d_l_p0=_b[hms_deep] - invttail(e(df_r),0.025)*_se[hms_deep]
		local d_u_p0=_b[hms_deep] + invttail(e(df_r),0.025)*_se[hms_deep]
		
		local d_b_p1=_b[L.hms_deep]
		local d_l_p1=_b[L.hms_deep] - invttail(e(df_r),0.025)*_se[L.hms_deep]
		local d_u_p1=_b[L.hms_deep] + invttail(e(df_r),0.025)*_se[L.hms_deep]
		
		local d_b_p2=_b[L2.hms_deep]
		local d_l_p2=_b[L2.hms_deep] - invttail(e(df_r),0.025)*_se[L2.hms_deep]
		local d_u_p2=_b[L2.hms_deep] + invttail(e(df_r),0.025)*_se[L2.hms_deep]
		
		* plot 
		preserve 
			clear 
			set obs 5
			
			gen eqtr=_n-3
			foreach var in d_b d_l d_u {
				gen `var'=. 
				replace `var' = ``var'_m2' if eqtr==-2 
				replace `var' = ``var'_m1' if eqtr==-1 
				replace `var' = ``var'_p0' if eqtr== 0 
				replace `var' = ``var'_p1' if eqtr== 1 
				replace `var' = ``var'_p2' if eqtr== 2 
			}
			
			foreach var in s_b s_l s_u {
				gen `var'=. 
				replace `var'=``var'' if eqtr==0
			}
			
			
			tw 	connected d_b eqtr, col(teal) lw(0.5) msymbol(circle_hollow) || ///
				line d_l d_u eqtr, col(gs9 gs9) lp(dash dash) || ///
				rcap s_l s_u eqtr, col(blue) || ///
				scatter s_b eqtr, col(blue) ///
				xtitle("Event year",size(huge)) ytitle("Employed per million",size(huge)) ///
				ylab(,labsize(vlarge) nogrid) yline(0,lp(dot) lw(0.6) lcol(black)) ///
				xlab(,labsize(vlarge)) ///
				legend(ring(0) pos(5) col(1) size(large) order(1 "dynamic" 5 "static")) ///
				graphregion(color(white)) xsize(4.5)
			gr export "$Rep_smokelabor/2_analysis/output_figures/figure3_c.pdf", replace 
		restore 
		
		
	}
	
	** panel d 
	if 1 {
		
		* static 
		reghdfe d_pmil_lau_lfp hms_deep [aw=seer_pop], a(fe_countyqtroy fe_styr) vce(cluster countyfip fe_stqtros)
		local s_b=_b[hms_deep]
		local s_l=_b[hms_deep] - invttail(e(df_r),0.025)*_se[hms_deep]
		local s_u=_b[hms_deep] + invttail(e(df_r),0.025)*_se[hms_deep]
		
		* dynamic 
		reghdfe d_pmil_lau_lfp F2.hms_deep F1.hms_deep hms_deep L1.hms_deep L2.hms_deep [aw=seer_pop], a(fe_countyqtroy fe_styr) vce(cluster countyfip fe_stqtros)
		local d_b_m2=_b[F2.hms_deep]
		local d_l_m2=_b[F2.hms_deep] - invttail(e(df_r),0.025)*_se[F2.hms_deep]
		local d_u_m2=_b[F2.hms_deep] + invttail(e(df_r),0.025)*_se[F2.hms_deep]
		
		local d_b_m1=_b[F.hms_deep]
		local d_l_m1=_b[F.hms_deep] - invttail(e(df_r),0.025)*_se[F.hms_deep]
		local d_u_m1=_b[F.hms_deep] + invttail(e(df_r),0.025)*_se[F.hms_deep]
		
		local d_b_p0=_b[hms_deep]
		local d_l_p0=_b[hms_deep] - invttail(e(df_r),0.025)*_se[hms_deep]
		local d_u_p0=_b[hms_deep] + invttail(e(df_r),0.025)*_se[hms_deep]
		
		local d_b_p1=_b[L.hms_deep]
		local d_l_p1=_b[L.hms_deep] - invttail(e(df_r),0.025)*_se[L.hms_deep]
		local d_u_p1=_b[L.hms_deep] + invttail(e(df_r),0.025)*_se[L.hms_deep]
		
		local d_b_p2=_b[L2.hms_deep]
		local d_l_p2=_b[L2.hms_deep] - invttail(e(df_r),0.025)*_se[L2.hms_deep]
		local d_u_p2=_b[L2.hms_deep] + invttail(e(df_r),0.025)*_se[L2.hms_deep]
		
		* plot 
		preserve 
			clear 
			set obs 5
			
			gen eqtr=_n-3
			foreach var in d_b d_l d_u {
				gen `var'=. 
				replace `var' = ``var'_m2' if eqtr==-2 
				replace `var' = ``var'_m1' if eqtr==-1 
				replace `var' = ``var'_p0' if eqtr== 0 
				replace `var' = ``var'_p1' if eqtr== 1 
				replace `var' = ``var'_p2' if eqtr== 2 
			}
			
			foreach var in s_b s_l s_u {
				gen `var'=. 
				replace `var'=``var'' if eqtr==0
			}
			
			
			tw 	connected d_b eqtr, col(teal) lw(0.5) msymbol(circle_hollow) || ///
				line d_l d_u eqtr, col(gs9 gs9) lp(dash dash) || ///
				rcap s_l s_u eqtr, col(blue) || ///
				scatter s_b eqtr, col(blue) ///
				xtitle("Event year",size(huge)) ytitle("LFP per million",size(huge)) ///
				ylab(,labsize(vlarge) nogrid) yline(0,lp(dot) lw(0.6) lcol(black)) ///
				xlab(,labsize(vlarge)) ///
				legend(ring(0) pos(5) col(1) size(large) order(1 "dynamic" 5 "static")) ///
				graphregion(color(white)) xsize(4.5)
			gr export "$Rep_smokelabor/2_analysis/output_figures/figure3_d.pdf", replace 
		restore 
	
	}
			
log close

